Parcial 2 TS YL
Yanina López
2023-11-12
#Librerias necesarias
library(astsa)
#install.packages(astsa)
#install.packages("ncdf4")
library(ncdf4)
#install.packages("data.table")
library(data.table)
#install.packages("forecast")
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo##
## Attaching package: 'forecast'## The following object is masked from 'package:astsa':
##
## gas#install.packages("urca")
library(urca)
#install.packages("rugarch")
library(rugarch)## Loading required package: parallel##
## Attaching package: 'rugarch'## The following object is masked from 'package:stats':
##
## sigma#install.packages("plotly")
library(plotly)## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutINTRODUCCIÓN
Modelar la temperatura del aire puede ofrecer una comprensión más profunda de los patrones climáticos, permitiéndo hacer predicciones más precisas y contribuir al entendimiento de los cambios climáticos a largo plazo y sus implicaciones en el ámbito de la meteorología, la climatología y otras disciplinas relacionadas.
El meridiano cero, al estar en Greenwich, longitud cero (polo norte) es un lugar de referencia para la hora estándar mundial (GMT/UTC). La recopilación de datos meteorológicos en este lugar proporciona información clave para la creación y mantenimiento de registros climatológicos globales, comprender los patrones climáticos globales y las variaciones estacionales. Estos datos son fundamentales para la investigación sobre el cambio climático y sus efectos.
La temperatura del aire en el meridiano cero es utilizada en modelos meteorológicos globales para prever patrones climáticos y eventos meteorológicos extremos en todo el mundo. Utilizar el meridiano cero como referencia para mediciones atmosféricas proporciona estándares internacionales consistentes para la comparación de datos meteorológicos en todo el mundo e incluso puede ser crucial para la aviación y la navegación, ya que afecta a la densidad del aire, lo que a su vez influye en la altitud de vuelo y en el rendimiento de las aeronaves.
La medición de la temperatura del aire a una atmósfera de presión en el meridiano cero también puede contribuir a estudios atmosféricos más amplios, como la comprensión de la estratificación térmica y la circulación atmosférica.
La National Oceanic and Atmospheric Administration (NOAA) posee satélites para capturar mediciones de variables utilizando la técnica de asimilación de datos que combina observaciones del mundo real con modelos numéricos para mejorar la precisión y la confiabilidad de las predicciones o estimaciones.
En ese sentido se obtuvo una serie de tiempo con datos mensuales de temperatura promedio del aire en escala Kelvin a 1000 milibares de presión que equivalen a aproximadamente una atmósfera de presión (sobre el nivel del mar) durante los últimos 44 años (Enero de 1979 - Diciembre de 2022).
FUENTE: (https://psl.noaa.gov/data/gridded/data.ncep.reanalysis2.html):
#cargue de datos
ruta_archivo <-"/Users/macbookair/Documents/DOCTORADO/SERIES DE TIEMPO/PARCIAL 2/air.mon.mean.nc"
Vwind <- nc_open(ruta_archivo)
#En caso de error con los datos descargue los datos desde el sgte enlace y reemplace la ruta de archivo: https://drive.google.com/file/d/1PGXFUVWz9pC9_rCmmmfuJgjxe0N9j65L/view?usp=sharing
attributes(Vwind$var)## $names
## [1] "time_bnds" "air"attributes(Vwind$dim)## $names
## [1] "lon" "lat" "level" "nbnds" "time"Se extraen los datos de temperatura del aire para lon:1,lat:1 y level:1: Temperatura promedio del aire en escala Kelvin a 1000 milibares de presión (que equivalen a aproximadamente una atmósfera de presión o sobre el nivel del mar), en la ubicación geografica de latitud cero sobre el meridiano de Greenwich.
mi_variable <-ncvar_get(Vwind, "air")
nc_close(Vwind)
Vwind.ts <- ts(mi_variable[1,1,1,c(1:537)], frequency = 12)
str(Vwind.ts)## Time-Series [1:537] from 1 to 45.7: 245 239 248 251 263 ...str(var)## function (x, y = NULL, na.rm = FALSE, use)Aplicando la metodología BOX-JENKINGS
IDENTIFICACIÓN: EDA
plot(decompose(Vwind.ts))
ggseasonplot(Vwind.ts)
ggseasonplot(Vwind.ts,polar=TRUE)
ggtsdisplay(Vwind.ts)
Modelo1<-sarima(Vwind.ts, p=1, d=1, q=1, P=1, D=1, Q=1, S=12)## initial value 1.363083
## iter 2 value 1.052293
## iter 3 value 0.992212
## iter 4 value 0.913991
## iter 5 value 0.903168
## iter 6 value 0.895102
## iter 7 value 0.893648
## iter 8 value 0.892104
## iter 9 value 0.892088
## iter 10 value 0.892035
## iter 11 value 0.892032
## iter 12 value 0.892032
## iter 13 value 0.892031
## iter 14 value 0.892031
## iter 15 value 0.892031
## iter 15 value 0.892031
## iter 15 value 0.892031
## final value 0.892031
## converged
## initial value 0.876949
## iter 2 value 0.869482
## iter 3 value 0.859587
## iter 4 value 0.858531
## iter 5 value 0.857976
## iter 6 value 0.857340
## iter 7 value 0.857290
## iter 8 value 0.857278
## iter 9 value 0.857277
## iter 10 value 0.857277
## iter 10 value 0.857277
## final value 0.857277
## converged
Modelo1$ttable## Estimate SE t.value p.value
## ar1 0.2040 0.0494 4.1340 0.0000
## ma1 -0.9687 0.0231 -41.9310 0.0000
## sar1 0.0620 0.0513 1.2077 0.2277
## sma1 -0.9140 0.0274 -33.4061 0.0000residuos<-Modelo1$fit$residuals
acf2(residuos)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF -0.01 0.06 -0.01 0.05 -0.03 0.03 -0.01 -0.02 -0.03 0.09 0.01 0.00 0.04
## PACF -0.01 0.06 -0.01 0.05 -0.03 0.03 -0.01 -0.03 -0.03 0.09 0.02 -0.01 0.04
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF 0.00 0.03 0.02 0.00 -0.03 0.02 0.03 -0.02 0 0.04 -0.01 0.01
## PACF -0.01 0.03 0.01 -0.01 -0.03 0.02 0.02 -0.02 0 0.03 -0.01 0.00
## [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
## ACF -0.01 -0.04 -0.08 -0.01 -0.02 0.00 -0.04 0.05 0.07 -0.05 0.07 -0.01
## PACF -0.02 -0.05 -0.07 -0.01 -0.02 0.01 -0.04 0.05 0.09 -0.07 0.06 0.00
## [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF -0.01 -0.02 0.02 -0.02 0.02 0.02 -0.06 -0.02 0.00 0.07 -0.02
## PACF -0.01 -0.01 0.02 0.00 0.03 0.03 -0.08 0.00 -0.01 0.08 -0.01Box.test(residuos, lag=8, type='Ljung-Box')##
## Box-Ljung test
##
## data: residuos
## X-squared = 4.9717, df = 8, p-value = 0.7606Box.test(residuos, lag=16, type='Ljung-Box')##
## Box-Ljung test
##
## data: residuos
## X-squared = 11.285, df = 16, p-value = 0.7915Box.test(residuos, lag=24, type='Ljung-Box')##
## Box-Ljung test
##
## data: residuos
## X-squared = 13.441, df = 24, p-value = 0.9583Debo verificar si hay raices unitarias en las estacionalidades y volatilidad. De la exploración de los datos, en las diferentes gráficas se puede identificar claramente un componente estacional en los datos originales por lo que una vez verificada la estacionariedad, el modelo propuesto debe poder representar la estacionalidad.
ESTIMACIÓN: Detección de raíces unitarias y verificar estabilidad
ESTACIONARIEDAD
Se verifica estacionariedad para serie original siguiendo el algoritmo de la prueba de Dickey Fuller aumentada:
ADF<-ur.df(Vwind.ts,type = "trend",lag=12)
summary(ADF)##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.3053 -1.2500 0.0843 1.4492 9.5314
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 93.202007 29.494369 3.160 0.00167 **
## z.lag.1 -0.363249 0.115071 -3.157 0.00169 **
## tt 0.003510 0.001379 2.546 0.01120 *
## z.diff.lag1 -0.300758 0.116838 -2.574 0.01033 *
## z.diff.lag2 -0.361666 0.111816 -3.234 0.00130 **
## z.diff.lag3 -0.470242 0.103546 -4.541 6.98e-06 ***
## z.diff.lag4 -0.455917 0.093954 -4.853 1.62e-06 ***
## z.diff.lag5 -0.528129 0.085185 -6.200 1.17e-09 ***
## z.diff.lag6 -0.527841 0.077591 -6.803 2.89e-11 ***
## z.diff.lag7 -0.557240 0.070112 -7.948 1.23e-14 ***
## z.diff.lag8 -0.619246 0.062223 -9.952 < 2e-16 ***
## z.diff.lag9 -0.692397 0.056998 -12.148 < 2e-16 ***
## z.diff.lag10 -0.622403 0.052706 -11.809 < 2e-16 ***
## z.diff.lag11 -0.422243 0.049733 -8.490 2.27e-16 ***
## z.diff.lag12 -0.080394 0.043949 -1.829 0.06795 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.501 on 509 degrees of freedom
## Multiple R-squared: 0.8399, Adjusted R-squared: 0.8355
## F-statistic: 190.8 on 14 and 509 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -3.1567 3.4028 5.0195
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.96 -3.41 -3.12
## phi2 6.09 4.68 4.03
## phi3 8.27 6.25 5.34Para evaluar la presencia de raíces unitarias en una serie temporal, lo que indica la no estacionariedad se seguirán los sigueintes pasos: Paso 1. Iniciar con el modelo menos restrictivo con un modelo que incluye tendencia,, usar \({\tau_3}\) para probar \({\gamma=0}\). Dado que la prueba tiene poca potencia si se rechaza entonces no se continua.
ADF1<-ur.df(Vwind.ts,type = "trend",lag=24, selectlags= "BIC")
summary(ADF1)##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.3048 -1.3445 0.1275 1.3024 9.4597
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 106.202015 29.879545 3.554 0.000415 ***
## z.lag.1 -0.414247 0.116625 -3.552 0.000419 ***
## tt 0.004223 0.001440 2.932 0.003519 **
## z.diff.lag1 -0.221842 0.114052 -1.945 0.052325 .
## z.diff.lag2 -0.271234 0.105311 -2.576 0.010296 *
## z.diff.lag3 -0.372543 0.094459 -3.944 9.16e-05 ***
## z.diff.lag4 -0.369191 0.085311 -4.328 1.82e-05 ***
## z.diff.lag5 -0.448534 0.076825 -5.838 9.51e-09 ***
## z.diff.lag6 -0.449656 0.068732 -6.542 1.51e-10 ***
## z.diff.lag7 -0.481309 0.059761 -8.054 5.96e-15 ***
## z.diff.lag8 -0.547546 0.052546 -10.420 < 2e-16 ***
## z.diff.lag9 -0.630295 0.045861 -13.743 < 2e-16 ***
## z.diff.lag10 -0.567666 0.044048 -12.887 < 2e-16 ***
## z.diff.lag11 -0.370945 0.041787 -8.877 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.523 on 498 degrees of freedom
## Multiple R-squared: 0.8374, Adjusted R-squared: 0.8331
## F-statistic: 197.3 on 13 and 498 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -3.552 4.2529 6.3239
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.96 -3.41 -3.12
## phi2 6.09 4.68 4.03
## phi3 8.27 6.25 5.34ADF2<-ur.df(Vwind.ts,type = "trend",lag=24, selectlags= "AIC")
summary(ADF2)##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.8099 -1.2605 0.0333 1.4008 8.9909
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 95.865549 34.377315 2.789 0.005501 **
## z.lag.1 -0.373916 0.134193 -2.786 0.005538 **
## tt 0.003898 0.001582 2.465 0.014056 *
## z.diff.lag1 -0.323768 0.135091 -2.397 0.016922 *
## z.diff.lag2 -0.330053 0.132635 -2.488 0.013164 *
## z.diff.lag3 -0.394377 0.130510 -3.022 0.002645 **
## z.diff.lag4 -0.361216 0.129402 -2.791 0.005454 **
## z.diff.lag5 -0.407340 0.128011 -3.182 0.001556 **
## z.diff.lag6 -0.386651 0.126576 -3.055 0.002377 **
## z.diff.lag7 -0.418297 0.125470 -3.334 0.000922 ***
## z.diff.lag8 -0.460856 0.124773 -3.694 0.000246 ***
## z.diff.lag9 -0.518100 0.124557 -4.160 3.77e-05 ***
## z.diff.lag10 -0.454651 0.124858 -3.641 0.000300 ***
## z.diff.lag11 -0.349479 0.125576 -2.783 0.005595 **
## z.diff.lag12 -0.110921 0.124937 -0.888 0.375080
## z.diff.lag13 -0.031793 0.120786 -0.263 0.792492
## z.diff.lag14 -0.136572 0.114844 -1.189 0.234943
## z.diff.lag15 -0.142753 0.107624 -1.326 0.185330
## z.diff.lag16 -0.166168 0.100754 -1.649 0.099742 .
## z.diff.lag17 -0.184907 0.093931 -1.969 0.049573 *
## z.diff.lag18 -0.232291 0.087293 -2.661 0.008048 **
## z.diff.lag19 -0.212242 0.079592 -2.667 0.007917 **
## z.diff.lag20 -0.196943 0.071869 -2.740 0.006364 **
## z.diff.lag21 -0.254496 0.062229 -4.090 5.06e-05 ***
## z.diff.lag22 -0.292448 0.054030 -5.413 9.77e-08 ***
## z.diff.lag23 -0.161028 0.044749 -3.598 0.000353 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.438 on 486 degrees of freedom
## Multiple R-squared: 0.8519, Adjusted R-squared: 0.8443
## F-statistic: 111.8 on 25 and 486 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -2.7864 2.7515 3.8821
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.96 -3.41 -3.12
## phi2 6.09 4.68 4.03
## phi3 8.27 6.25 5.34Teniendo en cuenta queValor del estadistico de prueba -2.7864 mayor que tau3 critico para cualquier nivel de significancia, por tanto no se rechaza Ho
Se selecciona el ajuste de ADF2 escogido por AIC con 23 rezagos o grados de libertad, cuyos residuos se comportan más como ruido blanco comparado con el ajuste ADF1 escogido por BIC con 11 rezagos. El summary de este ajuste arroja un valor del estadístico de -2.7864 el cual comparado con tau3 critico para cualquier nivel de significancia no permite rechazar la hipótesis nula. Este no es el resultado final ya que el proposito de esta parte es determinar el numero de lags. El pvalor dado en la tabla t para los coeficientes en particular el de gamma no es confiable porque la distribucion de los procesos con raiz unitaria no se comporta como t.
p-valor dado en la tabla t para los coeficientes, en particular el de gamma, no es confiable debido a que la distribución de los procesos con raíz unitaria no se comporta como una distribución t. La prueba de Dickey-Fuller es especialmente sensible a la cantidad de rezagos utilizados, y seleccionar el número óptimo de rezagos es una parte importante del proceso.
Se debe probar el modelo nuevamente según el paso 2 de la imagen para determinar si la serie no es estacionaria.
residuos<-ADF1@res
acf2(residuos)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF -0.03 0.02 0.04 0.05 0.03 0.03 0 0.00 0.01 0.01 -0.08 -0.18 0.00
## PACF -0.03 0.02 0.04 0.05 0.04 0.03 0 -0.01 0.00 0.01 -0.08 -0.19 -0.02
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF -0.06 0.05 0.06 0.03 -0.02 0.03 0.06 -0.04 -0.09 0.02 0.00 0.00
## PACF -0.05 0.07 0.10 0.07 0.00 0.03 0.05 -0.05 -0.12 -0.04 -0.05 -0.02
## [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## ACF -0.06 -0.06 -0.06 -0.01 -0.02 0.01 -0.01 0.05
## PACF -0.07 -0.01 -0.01 0.03 0.01 0.05 0.02 0.03residuos<-ADF2@res
acf2(residuos)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF 0 0.01 0.04 0.03 0.01 0.01 0.01 0 -0.01 -0.02 -0.02 -0.07 -0.06
## PACF 0 0.01 0.04 0.03 0.01 0.00 0.01 0 -0.01 -0.02 -0.03 -0.07 -0.06
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF 0.02 0.07 0.06 0.02 0.03 0.02 0.03 0.01 0.01 -0.08 -0.14 -0.04
## PACF 0.03 0.08 0.07 0.02 0.02 0.01 0.02 0.00 0.00 -0.10 -0.16 -0.05
## [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## ACF -0.06 -0.04 -0.05 -0.01 -0.01 0.04 -0.01 0.06
## PACF -0.05 -0.01 -0.02 0.01 0.00 0.05 0.01 0.07Se visualiza la autocorrelación de los residuos y evaluar su comportamiento.
# Hacemos la prueba de Dickey Fuller con el modelo restringido.
ADF_R<-ur.df(Vwind.ts,type = "drift",lag=23, selectlags= "Fixed")
summary(ADF_R)##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.747 -1.284 0.025 1.325 8.889
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.96201 16.08693 1.303 0.193173
## z.lag.1 -0.08061 0.06209 -1.298 0.194784
## z.diff.lag1 -0.60300 0.07376 -8.175 2.57e-15 ***
## z.diff.lag2 -0.59752 0.07653 -7.807 3.60e-14 ***
## z.diff.lag3 -0.64994 0.07937 -8.188 2.33e-15 ***
## z.diff.lag4 -0.60562 0.08337 -7.264 1.49e-12 ***
## z.diff.lag5 -0.64080 0.08644 -7.413 5.49e-13 ***
## z.diff.lag6 -0.60744 0.08977 -6.767 3.79e-11 ***
## z.diff.lag7 -0.62689 0.09297 -6.743 4.40e-11 ***
## z.diff.lag8 -0.65671 0.09651 -6.805 2.97e-11 ***
## z.diff.lag9 -0.70132 0.10034 -6.989 9.12e-12 ***
## z.diff.lag10 -0.62376 0.10465 -5.961 4.82e-09 ***
## z.diff.lag11 -0.50635 0.10859 -4.663 4.03e-06 ***
## z.diff.lag12 -0.25580 0.11032 -2.319 0.020822 *
## z.diff.lag13 -0.16559 0.10803 -1.533 0.125983
## z.diff.lag14 -0.25922 0.10368 -2.500 0.012742 *
## z.diff.lag15 -0.25192 0.09836 -2.561 0.010733 *
## z.diff.lag16 -0.26199 0.09316 -2.812 0.005118 **
## z.diff.lag17 -0.26794 0.08773 -3.054 0.002380 **
## z.diff.lag18 -0.30302 0.08249 -3.673 0.000266 ***
## z.diff.lag19 -0.27037 0.07605 -3.555 0.000414 ***
## z.diff.lag20 -0.24380 0.06927 -3.520 0.000472 ***
## z.diff.lag21 -0.29010 0.06063 -4.784 2.28e-06 ***
## z.diff.lag22 -0.31674 0.05303 -5.972 4.51e-09 ***
## z.diff.lag23 -0.17469 0.04459 -3.918 0.000102 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.448 on 488 degrees of freedom
## Multiple R-squared: 0.85, Adjusted R-squared: 0.8426
## F-statistic: 115.2 on 24 and 488 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -1.2983 1.0894
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1 6.43 4.59 3.78Paso 2. Verificar si debe incluirse la tendencia, usando \({\phi_3}\). Probemos entonces que \({a_2 = \gamma = 0}\). Recordemos las hipótesis involucradas y el cálculo de \({\phi_i}\)
Modelo restringido: \({H_0: \Delta Y_t = a_0+c_1 \Delta Y_{t-1}+...+c_p\Delta Y_{t-p} + e_t}\)
Modelo no restringido: \({H_1: \Delta Y_t = a_0+a_2t + \gamma \Delta Y_{t-1}+ c_1\Delta Y_{t-1}+...+c_p\Delta Y_{t-p} + e_t}\)
\({\phi_i = [(SSR_r - SSR_n)/r] / [SSR_n/(T-k)]}\)
\({SSR_r}\) = SSR restringido
\({SSR_n}\) = SSR no restringido
\({r}\) = Número de restricciones
\({T}\) = Número de observaciones utiles
\({k}\) = Número de parametros estimados en el modelo no restringido
\({T- k}\) = grados de libertad del modelo no restringido
SSR_U<-ADF2@testreg$sigma^2*ADF2@testreg$df[2]
SSR_R<-ADF_R@testreg$sigma^2*ADF_R@testreg$df[2]
r<-2 #Numero de restricciones
df_u<-ADF@testreg$df[2]
fi3=((SSR_R-SSR_U)/r)/(SSR_U/df_u)
fi3## [1] 3.182458Se calcula SSR para los modelos restringido y no restringido, define el número de restricciones y los grados de libertad, y luego calcula el estadístico \({\phi_3}\). con el fin de evaluar si la inclusión de la tendencia mejora significativamente el modelo, comparando el modelo no restringido (con tendencia) y el modelo restringido (sin tendencia).
Se espera que \({\phi_3}\) sea mayor a los puntos críticos del modelo no restringido para rechazar la hipótesis nula, se obtiene que es \({\phi_3}\) 3.038654 Al comparar con los puntos críticos del modelo ADF2 (NO RESTRINGIDO) se tiene que \({\phi_3}\) es menor por lo que por tanto no se rechaza H0, Quiere decir que entonces que no hay tendencia, puede haber raíz unitaria y debemos ir al siguiente paso del algoritmo.
# Miramos la ACF de los residuos para estar seguro que hay ruido blanco
residuos_ADF_R<-ADF_R@res
acf2(residuos_ADF_R)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF -0.01 0.01 0.03 0.03 0.01 0.01 0.01 0.01 -0.01 -0.02 -0.02 -0.06 -0.06
## PACF -0.01 0.01 0.03 0.03 0.01 0.00 0.01 0.01 -0.01 -0.02 -0.02 -0.07 -0.06
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF 0.02 0.07 0.05 0.01 0.02 0.01 0.03 0 0.01 -0.08 -0.15 -0.04
## PACF 0.02 0.07 0.06 0.02 0.02 0.01 0.02 0 0.00 -0.10 -0.17 -0.05
## [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## ACF -0.06 -0.04 -0.05 -0.01 -0.01 0.04 -0.01 0.06
## PACF -0.05 -0.01 -0.02 0.01 0.01 0.06 0.01 0.07Al revisar la ACF de los residuos vemos que hay una espiga que sobresale un poco y puede ser producto del azar. Por lo tanto, asumiremos que la ACF de los residuos es ruido blanco.
Paso 3. Estimaremos el modelo sin tendencia, es decir, \({\Delta Y_t = a_0 + \gamma * Y_{t-1}+c_1* \Delta Y{t-1} + ... + cp * \Delta Y_{t-p} + \epsilon_t}\). Para eso usamos \({\tau2}\) para verificar si hay raiz unitaria.
summary(ADF_R)##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.747 -1.284 0.025 1.325 8.889
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.96201 16.08693 1.303 0.193173
## z.lag.1 -0.08061 0.06209 -1.298 0.194784
## z.diff.lag1 -0.60300 0.07376 -8.175 2.57e-15 ***
## z.diff.lag2 -0.59752 0.07653 -7.807 3.60e-14 ***
## z.diff.lag3 -0.64994 0.07937 -8.188 2.33e-15 ***
## z.diff.lag4 -0.60562 0.08337 -7.264 1.49e-12 ***
## z.diff.lag5 -0.64080 0.08644 -7.413 5.49e-13 ***
## z.diff.lag6 -0.60744 0.08977 -6.767 3.79e-11 ***
## z.diff.lag7 -0.62689 0.09297 -6.743 4.40e-11 ***
## z.diff.lag8 -0.65671 0.09651 -6.805 2.97e-11 ***
## z.diff.lag9 -0.70132 0.10034 -6.989 9.12e-12 ***
## z.diff.lag10 -0.62376 0.10465 -5.961 4.82e-09 ***
## z.diff.lag11 -0.50635 0.10859 -4.663 4.03e-06 ***
## z.diff.lag12 -0.25580 0.11032 -2.319 0.020822 *
## z.diff.lag13 -0.16559 0.10803 -1.533 0.125983
## z.diff.lag14 -0.25922 0.10368 -2.500 0.012742 *
## z.diff.lag15 -0.25192 0.09836 -2.561 0.010733 *
## z.diff.lag16 -0.26199 0.09316 -2.812 0.005118 **
## z.diff.lag17 -0.26794 0.08773 -3.054 0.002380 **
## z.diff.lag18 -0.30302 0.08249 -3.673 0.000266 ***
## z.diff.lag19 -0.27037 0.07605 -3.555 0.000414 ***
## z.diff.lag20 -0.24380 0.06927 -3.520 0.000472 ***
## z.diff.lag21 -0.29010 0.06063 -4.784 2.28e-06 ***
## z.diff.lag22 -0.31674 0.05303 -5.972 4.51e-09 ***
## z.diff.lag23 -0.17469 0.04459 -3.918 0.000102 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.448 on 488 degrees of freedom
## Multiple R-squared: 0.85, Adjusted R-squared: 0.8426
## F-statistic: 115.2 on 24 and 488 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -1.2983 1.0894
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1 6.43 4.59 3.78Al comparar tau_2 se tiene que es -1,2983 con los valores críticos de tau2 y se tiene que patra para cualquier nivel de significancia es mayor por lo que no se rechaza la hipótesis nula. Lo que indica que aún es posible que exista una raíz unitaria. y debemos calcular phi1.
ADF_N<-ur.df(Vwind.ts,type = "none",lag=23, selectlags= "Fixed")
summary(ADF_N)##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.8253 -1.2291 0.0453 1.3413 8.6562
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 0.0002904 0.0004190 0.693 0.488653
## z.diff.lag1 -0.6796656 0.0445179 -15.267 < 2e-16 ***
## z.diff.lag2 -0.6706184 0.0520994 -12.872 < 2e-16 ***
## z.diff.lag3 -0.7194790 0.0587964 -12.237 < 2e-16 ***
## z.diff.lag4 -0.6716356 0.0662581 -10.137 < 2e-16 ***
## z.diff.lag5 -0.7034662 0.0718751 -9.787 < 2e-16 ***
## z.diff.lag6 -0.6667707 0.0774220 -8.612 < 2e-16 ***
## z.diff.lag7 -0.6832754 0.0823424 -8.298 1.04e-15 ***
## z.diff.lag8 -0.7101975 0.0874058 -8.125 3.68e-15 ***
## z.diff.lag9 -0.7519159 0.0925904 -8.121 3.80e-15 ***
## z.diff.lag10 -0.6712112 0.0981773 -6.837 2.42e-11 ***
## z.diff.lag11 -0.5514277 0.1030060 -5.353 1.33e-07 ***
## z.diff.lag12 -0.2982144 0.1054810 -2.827 0.004888 **
## z.diff.lag13 -0.2053551 0.1037045 -1.980 0.048242 *
## z.diff.lag14 -0.2961732 0.0997991 -2.968 0.003148 **
## z.diff.lag15 -0.2855764 0.0949771 -3.007 0.002776 **
## z.diff.lag16 -0.2923978 0.0902530 -3.240 0.001278 **
## z.diff.lag17 -0.2948816 0.0853173 -3.456 0.000595 ***
## z.diff.lag18 -0.3266302 0.0805314 -4.056 5.81e-05 ***
## z.diff.lag19 -0.2901430 0.0745697 -3.891 0.000114 ***
## z.diff.lag20 -0.2597082 0.0682293 -3.806 0.000159 ***
## z.diff.lag21 -0.3022332 0.0599564 -5.041 6.54e-07 ***
## z.diff.lag22 -0.3249980 0.0526919 -6.168 1.45e-09 ***
## z.diff.lag23 -0.1791509 0.0444872 -4.027 6.55e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.449 on 489 degrees of freedom
## Multiple R-squared: 0.8495, Adjusted R-squared: 0.8421
## F-statistic: 115 on 24 and 489 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: 0.693
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62SSR_U<-ADF_R@testreg$sigma^2*ADF_R@testreg$df[2]
SSR_R<-ADF_N@testreg$sigma^2*ADF_N@testreg$df[2]
r<-2 #Numero de restricciones
df_u<-ADF_R@testreg$df[2]
fi1=((SSR_R-SSR_U)/r)/(SSR_U/df_u)
fi1## [1] 0.8489648test_001<-fi1>ADF_R@cval[2,1]
test_005<-fi1>ADF_R@cval[2,2]
test_010<-fi1>ADF_R@cval[2,3]# Miramos la ACF de los residuos para estar seguro que hay ruido blanco
residuos_ADF_N<-ADF_N@res
acf2(residuos_ADF_N)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF -0.01 0.01 0.03 0.03 0.01 0 0.01 0.01 -0.01 -0.02 -0.02 -0.07 -0.06
## PACF -0.01 0.01 0.03 0.03 0.01 0 0.01 0.00 -0.01 -0.02 -0.02 -0.07 -0.06
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF 0.01 0.06 0.04 0.01 0.02 0.01 0.02 0.00 0.00 -0.09 -0.16 -0.04
## PACF 0.01 0.07 0.06 0.01 0.01 0.00 0.02 -0.01 -0.01 -0.10 -0.17 -0.06
## [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## ACF -0.06 -0.04 -0.06 -0.02 -0.01 0.04 -0.01 0.06
## PACF -0.06 -0.01 -0.03 0.00 0.00 0.05 0.00 0.07Teniendo en cuenta que phi_1 es menor para cualquier nivel de confianza que los puntos críticos del modelo no restringido para este punto, y los residuos se comportan como ruido blanco, entonces se acepta la hipótesis nula y se concluye que la serie de tiempo tiene tendencia estocástica por lo que hay raiz unitaria, sin intersecto ni tendencia.
DVwind<-(diff(Vwind.ts))
plot(DVwind)
Viendo el resultado que arroja la gráfica de la serie con una diferenciación, vemos un comportamiento que fluctua alrededor de 0.Ahora bien, para tener mayor certeza que la serie una vez diferenciada es estacionaria se aplica la función ndiffs que permite determinar el número de diferenciaciones que requiere la serie para volverse estacionaria, el resultado es 1:
ndiffs(Vwind.ts,alpha = 0.05,test = c("kpss", "adf", "pp"),type = c("level", "trend"),max.d = 2)## [1] 1Se revisan los residuos:
acf2(DVwind)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## ACF 0.55 0.28 -0.09 -0.35 -0.57 -0.61 -0.56 -0.37 -0.09 0.30 0.61 0.78
## PACF 0.55 -0.02 -0.34 -0.27 -0.32 -0.29 -0.35 -0.39 -0.44 -0.33 -0.17 0.16
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF 0.62 0.27 -0.07 -0.36 -0.54 -0.62 -0.54 -0.35 -0.08 0.27 0.62 0.76
## PACF 0.20 0.00 -0.03 -0.04 -0.01 -0.05 0.01 0.03 -0.06 -0.15 -0.03 0.09
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
## ACF 0.62 0.28 -0.07 -0.36 -0.53 -0.60 -0.53 -0.37 -0.06 0.28 0.58 0.76
## PACF 0.10 0.01 -0.04 -0.06 0.01 0.01 0.03 -0.03 0.00 -0.01 -0.10 0.09
## [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF 0.60 0.27 -0.08 -0.34 -0.53 -0.58 -0.52 -0.36 -0.07 0.26 0.61 0.71
## PACF 0.05 -0.04 -0.06 -0.01 -0.02 0.01 0.06 0.01 0.02 -0.09 0.06 0.00Para entender si hay patrones de correlación en los datos diferenciados se realiza el gráfico de autocorrelación de la serie diferenciada. En este se puede observar una autocorrelación Positiva en el Primer Rezago sugiere que hay una correlación positiva entre los valores actuales y los valores inmediatamente anteriores. En otras palabras, la temperatura del aire en un momento dado está correlacionada positivamente con la temperatura del aire en el momento anterior. Esto seguido de “ocho espigas hacia abajo” podrían indicar autocorrelaciones negativas en los rezagos siguientes. Esto podría interpretarse como una especie de patrón cíclico de cambios en la temperatura del aire. Ahora bien, los rezagos están por fuera de las bandas de confianza (las “bandas” generadas por la función acf()), podría indicar que algunas de estas autocorrelaciones son estadísticamente significativas.
Un patrón claro en las autocorrelaciones podría indica estacionalidad o repetición de patrones en la serie temporal de la temperatura del aire, lo que puede ser resultado de características específicas los datos de temperatura del aire, como factores climáticos, estacionales, o patrones diarios pueden influir en estos resultados. Se ratifica el comportamiento estacional de la serie aún despúes de una diferenciación y al no hallar ruido blanco implica que la diferenciación no hace que la serie sea estacionaria, en ese caso se se requiere probar una diferenciación estacional para cada 12 periodos:
plot(diff((Vwind.ts),lag=12))
ggseasonplot(diff(diff(Vwind.ts),lag=12))
ggseasonplot(diff(diff(Vwind.ts),lag=12), polar=TRUE)
ggsubseriesplot(diff(diff(Vwind.ts),lag=12))
acf2(diff(diff(Vwind.ts), lag=12))
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## ACF -0.43 -0.04 -0.06 0.07 -0.06 0.05 0.00 -0.05 -0.01 0.10 0.15 -0.46
## PACF -0.43 -0.28 -0.26 -0.14 -0.17 -0.10 -0.06 -0.12 -0.13 0.01 0.30 -0.29
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF 0.22 -0.03 0.06 0.00 0.02 -0.07 0.02 0.06 -0.03 -0.08 0.11 -0.07
## PACF -0.14 -0.17 -0.14 -0.04 -0.05 -0.06 -0.04 0.02 -0.03 -0.04 0.27 -0.16
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
## ACF 0.02 0.03 -0.01 -0.07 0.05 0.01 -0.01 -0.05 0.06 0.08 -0.16 0.14
## PACF -0.06 -0.06 -0.04 -0.04 -0.03 -0.05 -0.08 -0.06 -0.07 0.11 0.12 -0.03
## [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF -0.08 0.00 0.00 0.05 -0.08 0.06 0.04 -0.04 -0.04 -0.03 0.11 -0.11
## PACF -0.04 -0.09 -0.04 -0.05 -0.10 -0.06 0.04 -0.02 -0.04 0.02 0.14 -0.02Se hace la diferenciación estacional de 12 periodos al tratarse de datos mensuales
d12datos<-diff(Vwind.ts, lag=12)
acf2(d12datos)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF 0.18 0.05 0.01 0.05 -0.01 0.03 -0.01 -0.05 -0.01 0.04 -0.07 -0.43 -0.03
## PACF 0.18 0.02 -0.01 0.05 -0.03 0.03 -0.02 -0.05 0.02 0.04 -0.08 -0.42 0.14
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF -0.01 0.08 0.06 0.04 -0.01 0.04 0.05 -0.03 -0.07 0.01 -0.08 -0.04
## PACF 0.02 0.09 0.08 -0.02 0.01 0.02 0.01 -0.06 -0.01 0.00 -0.33 0.08
## [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
## ACF -0.03 -0.08 -0.11 -0.02 -0.01 -0.02 -0.01 0.08 0.07 -0.07 0.05 -0.05
## PACF -0.01 0.01 -0.01 0.00 -0.01 0.01 0.05 0.03 0.05 -0.13 -0.16 -0.02
## [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF -0.02 0.02 0.05 0.00 0.07 0.05 -0.04 -0.05 -0.01 0.08 0.00
## PACF -0.02 0.05 0.01 0.03 0.09 0.07 -0.04 0.02 0.05 -0.01 -0.13muestra<-window(d12datos, end=c(430,12))## Warning in window.default(x, ...): 'end' value not changedNótese que desaparecen las espigas estacionales. Se repite el proceso de estimación del número de rezagos en diferencia
TryMod<-ur.df(muestra, type="trend", lags=12, selectlags = "AIC")
summary(TryMod)##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.9712 -1.3831 0.0122 1.4974 9.9926
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.775e-02 2.464e-01 0.397 0.691798
## z.lag.1 -1.077e+00 1.287e-01 -8.368 5.96e-16 ***
## tt -5.024e-05 8.044e-04 -0.062 0.950225
## z.diff.lag1 2.674e-01 1.169e-01 2.287 0.022609 *
## z.diff.lag2 3.243e-01 1.129e-01 2.874 0.004230 **
## z.diff.lag3 2.995e-01 1.092e-01 2.743 0.006310 **
## z.diff.lag4 3.295e-01 1.042e-01 3.162 0.001660 **
## z.diff.lag5 2.995e-01 9.899e-02 3.026 0.002607 **
## z.diff.lag6 3.473e-01 9.350e-02 3.714 0.000227 ***
## z.diff.lag7 3.230e-01 8.808e-02 3.667 0.000272 ***
## z.diff.lag8 2.883e-01 8.096e-02 3.562 0.000404 ***
## z.diff.lag9 2.761e-01 7.361e-02 3.752 0.000196 ***
## z.diff.lag10 3.474e-01 6.467e-02 5.372 1.20e-07 ***
## z.diff.lag11 3.273e-01 5.628e-02 5.815 1.09e-08 ***
## z.diff.lag12 -1.273e-01 4.475e-02 -2.845 0.004625 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.687 on 497 degrees of freedom
## Multiple R-squared: 0.5424, Adjusted R-squared: 0.5295
## F-statistic: 42.08 on 14 and 497 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -8.3684 23.3628 35.0426
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.96 -3.41 -3.12
## phi2 6.09 4.68 4.03
## phi3 8.27 6.25 5.34resTryMod<-TryMod@res
acf2(resTryMod, max.lag = 30)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF 0 0 0.03 0.03 0 0 0.02 0.01 -0.03 -0.01 0.04 -0.17 -0.02
## PACF 0 0 0.03 0.03 0 0 0.02 0.01 -0.03 -0.01 0.04 -0.17 -0.01
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF 0.01 0.06 0.04 0 0.02 0.01 0.03 -0.04 0.03 0.00 -0.31 0
## PACF 0.01 0.07 0.05 0 0.01 0.00 0.04 -0.05 0.02 0.01 -0.35 0
## [,26] [,27] [,28] [,29] [,30]
## ACF -0.02 -0.03 -0.09 0 0.00
## PACF -0.02 0.01 -0.05 0 0.01Box.test(resTryMod, lag=8, type='Ljung-Box')##
## Box-Ljung test
##
## data: resTryMod
## X-squared = 1.3205, df = 8, p-value = 0.9953Box.test(resTryMod, lag=16, type='Ljung-Box')##
## Box-Ljung test
##
## data: resTryMod
## X-squared = 20.569, df = 16, p-value = 0.1957Box.test(resTryMod, lag=24, type='Ljung-Box')##
## Box-Ljung test
##
## data: resTryMod
## X-squared = 74.678, df = 24, p-value = 4.186e-07#Revisión de volatilidad
vol<-100*diff(log(Vwind.ts))
vol<-na.omit(vol)
plot(vol)
Una vez calculado y graficado los rendimientos porcentuales de la temperatura del aire, para obtener información sobre la volatilidad de la temperatura del aire a lo largo del tiempo. No se observan patrones notables en el gráfico, como picos o caídas bruscas, por lo que no existen indicios que puedan interpretarse como momentos de mayor o menor volatilidad en la temperatura del aire.
# Ajusta un modelo GARCH(1,1)
model_spec <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(0, 1)),
mean.model = list(armaOrder = c(0, 0)))
fit <- ugarchfit(spec = model_spec, data = d12datos)
# Obtiene los residuos estandarizados
residuals <- residuals(fit, standardize = TRUE)
# Crea un gráfico de dispersión condicional
plot(x = d12datos, y = residuals, type = "l",
xlab = "Tiempo", ylab = "Residuos Estandarizados",
main = "Gráfico de Dispersión Condicional")
# Añade una línea de referencia en y=0 para destacar los cambios en la volatilidad
abline(h = 0, col = "red", lty = 2)
Adicionalmente para analizar la volatilidad condicional en la serie de tiempo se realiza una gráfica de dispersión condicional que muestra la evolución de los residuos estandarizados a lo largo del tiempo que refleja que no hay mayor volatilidad debido a que gráficamente no se visualizan cambios destacados con respecto a la línea de referencia en y=0.
MODELACIÓN Y AJUSTE
A sabienda que existe estacionalidad en los datos originales, vamos a tratar de construir el modelo a partir de la misma, esto nos supone que en el modelo SARIMA \(d=1\).
Inicialmente se usa la función arima para ajustar un modelo ARIMA estacional. El orden del componente no estacional seleccionado como punto de partida según el análisis de autocorrelación de la serie de tiempo es (2,1,0) (indicando un AR(2) y una diferenciación de orden 1). El componente estacional se especifica como (0,0,5) (indicando un componente estacional de orden 5).
M1<-arima(muestra,order=c(2,1,0),seasonal=list(order=c(0,0,5)),fixed = c(0,NA,0,NA,0,0,NA))## Warning in arima(muestra, order = c(2, 1, 0), seasonal = list(order = c(0, :
## some AR parameters were fixed: setting transform.pars = FALSEM1##
## Call:
## arima(x = muestra, order = c(2, 1, 0), seasonal = list(order = c(0, 0, 5)),
## fixed = c(0, NA, 0, NA, 0, 0, NA))
##
## Coefficients:
## ar1 ar2 sma1 sma2 sma3 sma4 sma5
## 0 -0.0474 0 -0.1100 0 0 0.0537
## s.e. 0 0.0438 0 0.0493 0 0 0.0486
##
## sigma^2 estimated as 14.98: log likelihood = -1452.97, aic = 2913.94bic=AIC(M1,k = log(length(d12datos)))
bic## [1] 2930.994acf2(M1$residuals)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## ACF -0.44 0.00 -0.09 0.06 -0.05 0.05 -0.01 -0.04 0.01 0.08 0.17 -0.50
## PACF -0.44 -0.24 -0.25 -0.15 -0.17 -0.10 -0.07 -0.12 -0.11 0.02 0.34 -0.31
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF 0.23 -0.04 0.08 0.00 0.02 -0.06 0.03 0.05 -0.03 -0.08 0.07 0.01
## PACF -0.18 -0.17 -0.14 -0.04 -0.06 -0.06 -0.03 0.02 -0.02 -0.05 0.23 -0.14
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
## ACF 0.00 0.02 -0.02 -0.07 0.05 0.00 -0.01 -0.04 0.05 0.09 -0.14 0.09
## PACF -0.06 -0.05 -0.03 -0.04 -0.02 -0.05 -0.08 -0.04 -0.05 0.10 0.10 -0.02
## [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF -0.07 0.00 0.01 0.06 -0.08 0.06 0.04 -0.04 -0.02 -0.05 0.11 -0.09
## PACF -0.03 -0.07 -0.04 -0.04 -0.09 -0.06 0.05 -0.03 -0.03 0.03 0.10 -0.01acf2## function (series, max.lag = NULL, plot = TRUE, main = NULL, ylim = NULL,
## na.action = na.pass, ...)
## {
## acf = stats::acf
## pacf = stats::pacf
## frequency = stats::frequency
## num = length(series)
## xfreq = frequency(series)
## if (num > 59 & is.null(max.lag))
## max.lag = max(ceiling(10 + sqrt(num)), 4 * xfreq)
## if (num < 60 & is.null(max.lag))
## max.lag = floor(5 * log10(num + 5))
## if (max.lag > (num - 2))
## stop("Number of lags exceeds number of observations")
## if (is.null(main))
## main = paste("Series: ", deparse(substitute(series)))
## ACF = acf(series, max.lag, plot = FALSE, na.action = na.action,
## ...)$acf[-1]
## PACF = pacf(series, max.lag, plot = FALSE, na.action = na.action,
## ...)$acf
## LAG = (1:max.lag)/xfreq
## if (plot) {
## abline = graphics::abline
## par = graphics::par
## U = (-1/num) + (2/sqrt(num))
## L = (-1/num) - (2/sqrt(num))
## old.par <- par(no.readonly = TRUE)
## if (is.null(ylim)) {
## minA = min(ACF)
## maxA = max(ACF)
## minP = min(PACF)
## maxP = max(PACF)
## minu = min(minA, minP, L) - 0.01
## maxu = min(max(maxA + 0.1, maxP + 0.1), 1)
## ylim = c(minu, maxu)
## }
## par(mfrow = c(2, 1), cex.main = 1)
## Xlab = ifelse(xfreq > 1, paste("LAG ÷", xfreq), "LAG")
## tsplot(LAG, ACF, ylim = ylim, main = main, xlab = Xlab,
## ylab = "ACF", type = "h", ...)
## abline(h = c(0, L, U), lty = c(1, 2, 2), col = c(8, 4,
## 4))
## tsplot(LAG, PACF, ylim = ylim, main = NULL, xlab = Xlab,
## ylab = "PACF", type = "h", ...)
## abline(h = c(0, L, U), lty = c(1, 2, 2), col = c(8, 4,
## 4))
## on.exit(par(old.par))
## ACF <- round(ACF, 2)
## PACF <- round(PACF, 2)
## return(rbind(ACF, PACF))
## }
## else {
## return(cbind(ACF, PACF))
## }
## }
## <bytecode: 0x7fbdc75b3d78>
## <environment: namespace:astsa>En este primer modelo M1 se obtiene que no hay ruido blanco por lo que se deben explorar otras opciones para ajustar el modelo.
Se usa la función autoarima para obtener aproximación de parámetros en la modelación:
auto.arima(muestra)## Series: muestra
## ARIMA(2,0,1)(0,0,1)[12] with non-zero mean
##
## Coefficients:
## ar1 ar2 ma1 sma1 mean
## -0.5299 0.2388 0.7587 -0.8854 0.1154
## s.e. 0.1813 0.0498 0.1815 0.0291 0.0195
##
## sigma^2 = 5.337: log likelihood = -1191.2
## AIC=2394.4 AICc=2394.56 BIC=2419.98M2<-arima(muestra,order=c(2,0,1),seasonal=list(order=c(0,1,1),period=12))
M2##
## Call:
## arima(x = muestra, order = c(2, 0, 1), seasonal = list(order = c(0, 1, 1), period = 12))
##
## Coefficients:
## ar1 ar2 ma1 sma1
## 0.6416 -0.0587 -0.4695 -1.0000
## s.e. 0.7543 0.1486 0.7525 0.0172
##
## sigma^2 estimated as 9.132: log likelihood = -1317.93, aic = 2645.85bic=AIC(M2,k = log(length(muestra)))
bic## [1] 2667.17acf2(M2$residuals)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF 0 0 -0.02 0.05 -0.03 0.03 -0.01 -0.05 0 0.07 0 -0.43 0.04
## PACF 0 0 -0.02 0.05 -0.03 0.03 -0.01 -0.05 0 0.07 0 -0.43 0.06
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF -0.01 0.08 0.05 0.03 -0.03 0.03 0.06 -0.03 -0.07 0.04 -0.06 -0.02
## PACF 0.01 0.08 0.10 -0.01 0.00 0.01 0.03 -0.04 -0.02 0.06 -0.31 0.04
## [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
## ACF -0.01 -0.05 -0.09 0.00 -0.01 -0.01 -0.01 0.08 0.07 -0.09 0.08 -0.06
## PACF -0.02 0.00 -0.01 -0.01 -0.02 -0.01 0.06 0.04 0.08 -0.07 -0.13 -0.04
## [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF -0.02 0.01 0.05 -0.02 0.07 0.05 -0.04 -0.05 -0.01 0.09 -0.03
## PACF -0.04 0.03 0.02 0.02 0.06 0.06 -0.01 0.02 0.06 0.03 -0.07acf2## function (series, max.lag = NULL, plot = TRUE, main = NULL, ylim = NULL,
## na.action = na.pass, ...)
## {
## acf = stats::acf
## pacf = stats::pacf
## frequency = stats::frequency
## num = length(series)
## xfreq = frequency(series)
## if (num > 59 & is.null(max.lag))
## max.lag = max(ceiling(10 + sqrt(num)), 4 * xfreq)
## if (num < 60 & is.null(max.lag))
## max.lag = floor(5 * log10(num + 5))
## if (max.lag > (num - 2))
## stop("Number of lags exceeds number of observations")
## if (is.null(main))
## main = paste("Series: ", deparse(substitute(series)))
## ACF = acf(series, max.lag, plot = FALSE, na.action = na.action,
## ...)$acf[-1]
## PACF = pacf(series, max.lag, plot = FALSE, na.action = na.action,
## ...)$acf
## LAG = (1:max.lag)/xfreq
## if (plot) {
## abline = graphics::abline
## par = graphics::par
## U = (-1/num) + (2/sqrt(num))
## L = (-1/num) - (2/sqrt(num))
## old.par <- par(no.readonly = TRUE)
## if (is.null(ylim)) {
## minA = min(ACF)
## maxA = max(ACF)
## minP = min(PACF)
## maxP = max(PACF)
## minu = min(minA, minP, L) - 0.01
## maxu = min(max(maxA + 0.1, maxP + 0.1), 1)
## ylim = c(minu, maxu)
## }
## par(mfrow = c(2, 1), cex.main = 1)
## Xlab = ifelse(xfreq > 1, paste("LAG ÷", xfreq), "LAG")
## tsplot(LAG, ACF, ylim = ylim, main = main, xlab = Xlab,
## ylab = "ACF", type = "h", ...)
## abline(h = c(0, L, U), lty = c(1, 2, 2), col = c(8, 4,
## 4))
## tsplot(LAG, PACF, ylim = ylim, main = NULL, xlab = Xlab,
## ylab = "PACF", type = "h", ...)
## abline(h = c(0, L, U), lty = c(1, 2, 2), col = c(8, 4,
## 4))
## on.exit(par(old.par))
## ACF <- round(ACF, 2)
## PACF <- round(PACF, 2)
## return(rbind(ACF, PACF))
## }
## else {
## return(cbind(ACF, PACF))
## }
## }
## <bytecode: 0x7fbdc75b3d78>
## <environment: namespace:astsa>M3<-sarima(muestra,p=2,d=0,q=1,P=0,D=1,Q=4,S=12)## initial value 1.636228
## iter 2 value 1.300917
## iter 3 value 1.138853
## iter 4 value 1.118309
## iter 5 value 1.102626
## iter 6 value 1.081712
## iter 7 value 1.074912
## iter 8 value 1.068452
## iter 9 value 1.065561
## iter 10 value 1.064848
## iter 11 value 1.063966
## iter 12 value 1.063702
## iter 13 value 1.063667
## iter 14 value 1.063600
## iter 15 value 1.063399
## iter 16 value 1.062820
## iter 17 value 1.062453
## iter 18 value 1.061801
## iter 19 value 1.061584
## iter 20 value 1.061167
## iter 21 value 1.061067
## iter 22 value 1.060661
## iter 23 value 1.060520
## iter 24 value 1.060175
## iter 25 value 1.060071
## iter 26 value 1.060058
## iter 27 value 1.060014
## iter 28 value 1.059982
## iter 29 value 1.059933
## iter 30 value 1.059915
## iter 31 value 1.059911
## iter 32 value 1.059911
## iter 32 value 1.059911
## iter 32 value 1.059911
## final value 1.059911
## converged
## initial value 1.028001
## iter 2 value 1.009554
## iter 3 value 1.006791
## iter 4 value 0.986489
## iter 5 value 0.981427
## iter 6 value 0.949323
## iter 7 value 0.946465
## iter 8 value 0.945083
## iter 9 value 0.940426
## iter 10 value 0.939957
## iter 11 value 0.939227
## iter 12 value 0.938873
## iter 13 value 0.938389
## iter 14 value 0.938212
## iter 15 value 0.938199
## iter 16 value 0.938198
## iter 17 value 0.938197
## iter 18 value 0.938194
## iter 18 value 0.938194
## final value 0.938194
## converged
M3$ttable## Estimate SE t.value p.value
## ar1 0.6895 1.1201 0.6156 0.5384
## ar2 -0.0455 0.2721 -0.1671 0.8674
## ma1 -0.4902 1.1208 -0.4374 0.6620
## sma1 -1.8247 0.1530 -11.9232 0.0000
## sma2 0.7507 0.1551 4.8415 0.0000
## sma3 0.0899 0.0998 0.9010 0.3680
## sma4 -0.0137 0.0488 -0.2817 0.7783
## constant 0.0004 0.0002 2.1393 0.0329M3$AIC## [1] 4.749353M3$BIC## [1] 4.823744Para verificar si los residuos del modelo son ruido blanco, es decir validar si no hay patrones de autocorrelación se aplica la prueba Ljung-Box, con el fin de indentificar si la serie de tiempo exhibe autocorrelación significativa en varios lags. En este caso, los términos sma1 y sma2 son significativos, lo que sugiere que tienen un impacto estadísticamente significativo en el modelo.
Se tiene que el modelo M3 presenta buen ajuste, los coeficientes del modelo según la tabla t son significativos y los residuos del mismo se comportan como ruido blanco y presenta el mejor AIC.
CHEQUEO Y SELECCIÓN DEL MODELO:
Teniendo en cuenta que La hipótesis nula (H0) en la prueba de Box-Ljung es que no hay autocorrelación significativa en los residuos de la serie de tiempo hasta un cierto número de lags y La hipótesis alternativa (H1) es que hay autocorrelación significativa en los residuos hasta un cierto número de lags. La prueba calcula una estadística de prueba, conocida como estadístico Q de Ljung-Box, que se basa en los valores autocorrelacionados de los residuos hasta un número específico de lags. Se tiene que el modelo M4 presenta buen ajuste para la serie de tiempo ya que o hay autocorrelación significativa en los residuos de la serie de tiempo hasta un cierto número de lags y estos son ruido blanco.
Ahora bien alcomparar el BIC (Criterio de Información Bayesiano) y el AIC (Criterio de Información de Akaike) favorecen modelos más pequeños y parsimoniosos pero el BIC tiende a penalizar la complejidad más fuertemente que el AIC de los dos modelos se selecciona el modelo M3 por tener mejor BIC ``
Se puede apreciar gráficamente que el modelo M3 presenta un comportamiento similar a los datos de temperatura del aire, es decir que representa la serie de tiempo de temperatura media del aire de los últimos 44 años a una atmósfera de presión en el meridiando Greenwich, latitud cero.
Se aprecia estacionalidad año tras año si mayores fluctuaciones en los rendimientos año tras año teniendo en cuenta que los comportamientos mensuales a través del tiempo son bastante estables por lo tanto no se contemplará modelar la volatilidad
M4<-sarima(muestra, p=2,d=0,q=1,P=0,D=1,Q=4,S=12)## initial value 1.636228
## iter 2 value 1.300917
## iter 3 value 1.138853
## iter 4 value 1.118309
## iter 5 value 1.102626
## iter 6 value 1.081712
## iter 7 value 1.074912
## iter 8 value 1.068452
## iter 9 value 1.065561
## iter 10 value 1.064848
## iter 11 value 1.063966
## iter 12 value 1.063702
## iter 13 value 1.063667
## iter 14 value 1.063600
## iter 15 value 1.063399
## iter 16 value 1.062820
## iter 17 value 1.062453
## iter 18 value 1.061801
## iter 19 value 1.061584
## iter 20 value 1.061167
## iter 21 value 1.061067
## iter 22 value 1.060661
## iter 23 value 1.060520
## iter 24 value 1.060175
## iter 25 value 1.060071
## iter 26 value 1.060058
## iter 27 value 1.060014
## iter 28 value 1.059982
## iter 29 value 1.059933
## iter 30 value 1.059915
## iter 31 value 1.059911
## iter 32 value 1.059911
## iter 32 value 1.059911
## iter 32 value 1.059911
## final value 1.059911
## converged
## initial value 1.028001
## iter 2 value 1.009554
## iter 3 value 1.006791
## iter 4 value 0.986489
## iter 5 value 0.981427
## iter 6 value 0.949323
## iter 7 value 0.946465
## iter 8 value 0.945083
## iter 9 value 0.940426
## iter 10 value 0.939957
## iter 11 value 0.939227
## iter 12 value 0.938873
## iter 13 value 0.938389
## iter 14 value 0.938212
## iter 15 value 0.938199
## iter 16 value 0.938198
## iter 17 value 0.938197
## iter 18 value 0.938194
## iter 18 value 0.938194
## final value 0.938194
## converged
M4$ttable## Estimate SE t.value p.value
## ar1 0.6895 1.1201 0.6156 0.5384
## ar2 -0.0455 0.2721 -0.1671 0.8674
## ma1 -0.4902 1.1208 -0.4374 0.6620
## sma1 -1.8247 0.1530 -11.9232 0.0000
## sma2 0.7507 0.1551 4.8415 0.0000
## sma3 0.0899 0.0998 0.9010 0.3680
## sma4 -0.0137 0.0488 -0.2817 0.7783
## constant 0.0004 0.0002 2.1393 0.0329residuos4<-M4$fit$residuals
acf2(residuos4)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF 0 0.01 -0.03 0.04 -0.03 0.04 -0.01 -0.02 -0.02 0.09 0.01 -0.01 0.03
## PACF 0 0.01 -0.03 0.04 -0.03 0.04 0.00 -0.02 -0.01 0.08 0.01 -0.01 0.04
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF -0.01 0.04 0.03 0 -0.03 0.02 0.04 -0.02 -0.01 0.03 0.01 0
## PACF -0.01 0.05 0.03 0 -0.02 0.02 0.04 -0.03 0.00 0.02 0.01 0
## [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
## ACF -0.01 -0.05 -0.07 0.00 -0.01 0 -0.02 0.06 0.07 -0.06 0.06 -0.01
## PACF -0.02 -0.05 -0.06 -0.01 -0.02 0 -0.02 0.05 0.08 -0.07 0.07 0.00
## [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF -0.02 -0.02 0.03 -0.01 0.02 0.02 -0.06 -0.03 0.00 0.06 -0.02
## PACF -0.02 -0.02 0.02 0.00 0.03 0.02 -0.08 -0.01 -0.01 0.05 -0.01Box.test(residuos4, lag=8, type='Ljung-Box')##
## Box-Ljung test
##
## data: residuos4
## X-squared = 3.0098, df = 8, p-value = 0.9337Box.test(residuos4, lag=16, type='Ljung-Box')##
## Box-Ljung test
##
## data: residuos4
## X-squared = 9.5372, df = 16, p-value = 0.8897Box.test(residuos4, lag=24, type='Ljung-Box')##
## Box-Ljung test
##
## data: residuos4
## X-squared = 11.871, df = 24, p-value = 0.9813según estos resultados del test de Ljung-Box, no hay evidencia significativa de autocorrelación en los residuos del modelo SARIMA hasta los rezagos considerados (8, 16, y 24). Esto sugiere que el modelo SARIMA ha capturado adecuadamente la estructura de autocorrelación en la serie de tiempo. Lo anterior implica que se ha encontrado el Modelo.
Pronóstico
A continuación se presenta el resultado del análisis de la serie de tiempo de temperatura del aire, se identificó que la serie tiene un comportamiento estacional inicialmente no estacionaria que tuvo que diferenciarse estacionalmente para poder ser modelada, no fue necesario modelar volatilidad ni cambios estructurales en la serie y con la modelación se pudo hallar un modelo SARIMA con un componente autoregresivo (AR), un componente de media móvil (MA), tanto en la parte no estacional como en la estacional, y también incluye diferenciación para manejar tendencias no estacionarias y estacionalidad mensual. La elección específica de estos parámetros se basa en el análisis y la comprensión de la serie de tiempo particular que se ha modelado.
A continuación se presenta la gráfica del modelo con los datos de muestra frente a la serie original de temperatura del aire.
sarima.for(muestra,n.ahead=24,p=2,d=0,q=1,P=0,D=1,Q=4,S=12)## $pred
## Jan Feb Mar Apr May Jun Jul
## 45
## 46 0.7767042 4.1919155 3.1136137 5.9711904 3.6871036 2.0276635 1.8500717
## 47 0.3802075 0.7470753 0.6233433 0.7555095 0.5028296 0.3187339 0.2930456
## Aug Sep Oct Nov Dec
## 45 1.4909667 3.6729078 2.0826290
## 46 1.2750707 0.5939942 0.6005659 0.7797546 0.5516416
## 47 0.2752712 0.2099281
##
## $se
## Jan Feb Mar Apr May Jun Jul Aug
## 45
## 46 2.387892 2.389150 2.389627 2.389807 2.389876 2.389901 2.389911 2.389914
## 47 3.065040 3.065670 3.065909 3.066000 3.066034 3.066047 3.066052 3.066054
## Sep Oct Nov Dec
## 45 2.330450 2.376299 2.385897
## 46 2.389914 3.035797 3.058569 3.063387
## 47 3.066055lines(Vwind.ts)
fig <- plot_ly(y = Vwind.ts, name = "Original Data", type = 'scatter', mode = 'lines')
fig <- fig %>% add_trace(y = Vwind.ts+M4$fit$residuals, name = "Predicted Values", connectgaps = TRUE)
figSe hace la comparación final de los datos originales vs los datos generados por el modelo.